!*******************************************************************************
!This module calculates the value functions in tenure 8 by value function 
!iteration and the value functions in tenures 1 to 7 by backward induction
!NOTE: estimates in the paper are not reported scaled by (1-delta)
!*******************************************************************************
MODULE mod_values


USE mod_grid
USE mod_param_initial
USE mod_parameters
USE mod_update
USE mod_types


IMPLICIT NONE


CONTAINS

!***************************************************
!(A) Program to compute values by value iteration
!***************************************************
! Structure for value iteration:
!    * value = PBAR
!    * it = 0  
!    * diff = constv
!    * DO 
!       value_old = value
!       ...
!       value = f(value_old)
!       ...
!       diff = MAXVAL(DABS(value-value_old))
!       it = it+1
!       IF (diff<=tol_level .or. it>max_iter) EXIT  
!      END DO


FUNCTION argu3(task,phi,svalue) 


REAL(dp), INTENT(IN):: phi, svalue
INTEGER, INTENT(IN)::  task
REAL(dp), PARAMETER :: vvtinys = 1.0d-300, vvtiny2 = 1.0d-12, & 
                       maxmax= 0.0_dp, cconst = 1.0_dp !valmax= 1300.0_dp
REAL(dp):: argu3




IF (task==1) THEN

argu3 = -maxmax+ svalue


ELSE IF (task==2) THEN

argu3 = -maxmax+ svalue



ELSE IF (task>=3) THEN


argu3 = -maxmax+ svalue



END IF




END FUNCTION argu3




FUNCTION drift(drift_level,drift_tenure,t2_L1,t2_L2, t3_L1, t3_L2, t3_L3, t4_L1, t4_L2, t4_L3, t5_L1, t5_L2, t5_L3,&
               t6_L1, t6_L2, t6_L3)


INTEGER, INTENT(IN):: drift_level,drift_tenure, t2_L1, t2_L2, t3_L1, t3_L2, t3_L3, t4_L1, t4_L2, t4_L3,t5_L1, t5_L2, t5_L3, &
                      t6_L1, t6_L2, t6_L3

REAL(dp):: drift


IF(drift_level==1) THEN



IF(drift_tenure==1) THEN

  drift = 1000.0_dp


 ELSE IF(drift_tenure==2) THEN
  
  
  drift =  0.0_dp



  ELSE IF(drift_tenure==3) THEN
    
 
  drift = 0.0_dp   


  ELSE IF(drift_tenure==4) THEN
   

 
  drift = 0.0_dp -d*( REAL(t3_L2,dp) )

  
  ELSE IF(drift_tenure==5) THEN
  
 
 
  drift = 0.0_dp -d32*( REAL(t4_L2,dp) )



  ELSE IF(drift_tenure==6) THEN
  
  
  
  drift = 0.0_dp -d32*( REAL(t5_L2,dp) )
   
        


  ELSE IF(drift_tenure==7) THEN
  
 
  drift = 0.0_dp  -d32*( REAL(t6_L2,dp) )  
    

 


  ELSE IF(drift_tenure==8) THEN
  
  
  drift = 0.0_dp 



 
 END IF



ELSE IF (drift_level==2) THEN


 IF(drift_tenure==1) THEN

  drift = 0.0_dp


  ELSE IF(drift_tenure==2) THEN
 
  
  drift =   0.0_dp   


  ELSE IF(drift_tenure==3) THEN
  
  
  drift =  0.0_dp 


  ELSE IF(drift_tenure==4) THEN
  
  
  drift =   0.0_dp 


  ELSE IF(drift_tenure==5) THEN
  
    
  drift =  0.0_dp 
        

  ELSE IF(drift_tenure==6) THEN
  

   drift = 0.0_dp 
         

  
   ELSE IF(drift_tenure==7) THEN
  
  
  drift =  0.0_dp 
         

   ELSE IF(drift_tenure==8) THEN
  
 
  
  drift = 0.0_dp
 
 END IF




ELSE IF (drift_level>=3) THEN


IF(drift_tenure==1) THEN
  
  
  drift = 0.0_dp


  ELSE IF(drift_tenure==2) THEN
  

  drift = 0.0_dp


  ELSE IF(drift_tenure==3) THEN
  
  
  drift =  0.0_dp


  ELSE IF(drift_tenure==4) THEN
  
  
  drift = d31*( REAL(t3_L3,dp) )
 
  
  
  ELSE IF(drift_tenure==5) THEN
  
   drift =  d30*( REAL(t4_L3,dp) )  
    
 
  ELSE IF(drift_tenure==6) THEN
  
  drift = d22*( REAL(t5_L3,dp) )  
  
 
  ELSE IF(drift_tenure==7) THEN

  
  drift =  d24*( REAL(t6_L3,dp) ) 
  
 
  ELSE IF(drift_tenure==8) THEN
  
    
   drift = 0.0_dp 

 
 END IF

END IF



END FUNCTION drift



SUBROUTINE value_infinite(value8,task,V_0,V_11, V_21, V_31, V_12, V_22, V_32, V_13, V_23, V_33, &
     V_14, V_24, V_34, V_15, V_25, V_35, V_16, V_26, V_36, V_17, V_27, V_37, V_18, V_28, V_38, phi_1_grid_point,phi_1_grid_point_max, phi_1_star,phi_2_grid_point,&
	 phi_2_star,phi_3_star)

IMPLICIT NONE


 
INTEGER, DIMENSION (phi_grid_size,3,2):: min_phi

INTEGER:: tenure3_L1, tenure3_L2, tenure3_L3, tenure4_L1, tenure4_L2, tenure4_L3, tenure5_L1, tenure5_L2, tenure5_L3, tenure6_L1, tenure6_L2, tenure6_L3
INTEGER:: vl3, vl4,vl5,vl6
INTEGER:: vi1, vvi1, vi2, vi3, it, phi_1_grid_point, phi_1_grid_point_max, phi_2_grid_point
INTEGER:: min_omega1high(1), min_omega2high(1), min_omega3high(1), min_omega1low(1), min_omega2low(1), min_omega3low(1)
INTEGER:: min_phi_1_star(1), max_phi_1_star(1), &
min_phi_2_star(1), min_phi_3_star(1), max_task(1)
REAL(dp), DIMENSION (phi_grid_size):: value, value_old, value0, value0_old, value1_old, value2_old, value3_old, &
    diff_phi_12, diff_phi_11, diff_phi_22, diff_phi_21, diff_phi_32, diff_phi_31, diff_omega1high, diff_omega2high, &
	diff_omega3high, diff_omega1low, diff_omega2low, diff_omega3low
	
	
REAL(dp), DIMENSION(phi_grid_size):: V_18, V_28, V_38, V_8, value8, value8_old, V_18_old, V_28_old, V_38_old

REAL(dp), DIMENSION(phi_grid_size ,3,3,3,3):: V_17, V_27, V_37, V_7, value7

REAL(dp), DIMENSION(phi_grid_size ,3,3,3):: V_16, V_26, V_36, V_6, value6

REAL(dp), DIMENSION(phi_grid_size ,3,3):: V_15, V_25, V_35, V_5, value5

REAL(dp), DIMENSION(phi_grid_size ,3):: V_14, V_24, V_34, V_4, value4

REAL(dp), DIMENSION(phi_grid_size ):: V_13, V_23, V_33, V_3, value3

REAL(dp), DIMENSION(phi_grid_size ):: V_12, V_22, V_32, V_2, value2

REAL(dp), DIMENSION(phi_grid_size ):: V_11, V_21, V_31, V_1, value1

REAL(dp), DIMENSION (phi_grid_size)::  V_0, V_0_old


INTEGER, DIMENSION (phi_grid_size):: task


INTEGER, PARAMETER:: max_iter = 3000


!REAL(dp), PARAMETER:: tol_level= 30
REAL(dp), PARAMETER:: tol_level=1.0d-04, correct5 = 1000.0_dp, correct4 = 1000.0_dp, correct3 = 1000.0_dp, &
                      correct2 = 1000.0_dp, correct1 = 1000.0_dp, correct = 500.0_dp
                      
REAL(dp):: trydiff

REAL(dp):: diff,  diff0, diff1, diff2, diff3, max_diff, value_new, value0_new, value1_new, value2_new, value3_new, &
V_0_operator, V_1_operator, V_2_operator, &
V_3_operator, phi_1_star, phi_2_star, phi_2_star_temp, phi_3_star, &
V_0_new, V_1_new, V_2_new, V_3_new, V_00_new, V_10_new, V_20_new, V_30_new, V_01_new, V_11_new, V_21_new, V_31_new, &
V_02_new, V_12_new, V_22_new, V_32_new, V_03_new, V_13_new, V_23_new, V_33_new, &
value8_new, V_18_new, V_28_new, V_38_new, try
REAL(dp), DIMENSION(4):: B_operator
REAL(dp), DIMENSION(4):: B0_operator, B1_operator, B2_operator, B3_operator

REAL(dp), PARAMETER:: vtiny = 1.0d-300, vtiny2 = 1.0d-20 

INTEGER, DIMENSION(:), ALLOCATABLE:: ind_task_1
INTEGER:: phi_grid_ind
!REAL(dp), PARAMETER:: max_exp = 550.00_dp
REAL(dp), PARAMETER:: max_exp = 0.00_dp, pen= 0.0_dp

REAL(dp), PARAMETER:: constv = 1.00_dp
REAL(dp):: arg3(3)
!REAL(dp), PARAMETER:: tau=1.0_dp


REAL(dp), PARAMETER:: largeno=1000.0_dp



CALL grid_generation  


!PRINT*, 'age1', age1

DO vi2 = 1, phi_grid_size  
   !diff_omega1high(vi2) = DABS(omega_11-phi_grid(vi2))
   !diff_omega2high(vi2) = DABS(omega_12-phi_grid(vi2))
   !diff_omega3high(vi2) = DABS(omega_13-phi_grid(vi2))
   diff_omega1low(vi2) = DABS(omega_21 - phi_grid(vi2))
   diff_omega2low(vi2) = DABS(omega_22 - phi_grid(vi2))
   diff_omega3low(vi2) = DABS(omega_23 - phi_grid(vi2))
END DO

!min_omega1high = MINLOC(diff_omega1high)  
!min_omega2high = MINLOC(diff_omega2high)  
!min_omega3high = MINLOC(diff_omega3high)  

min_omega1low = MINLOC(diff_omega1low)  
min_omega2low = MINLOC(diff_omega2low)  
min_omega3low = MINLOC(diff_omega3low)  


value8 = vtiny2

  V_0  = vtiny2
  
  V_18 = vtiny2
 

  V_28 = vtiny2


  V_38 = vtiny2
  


CALL grid_generation  


CALL sub_min_phi(min_phi)

!***************************************************************
!(1) Initialize iteration index, value difference and value
!***************************************************************
it = 0           
diff = 10000_dp



!********************
!(2) Compute values
!********************
!--------------------------
!(a) Loop for grid search
!--------------------------
DO 

value8_old = value8

V_0_old = V_0

V_18_old = V_18

V_28_old = V_28

V_38_old = V_38



!-------------------------------------------
!(b) Loop for the current period prior
!-------------------------------------------
DO vi1 = 1, phi_grid_size 


!********************************
!(3) Compute the value function 
!********************************
V_0_new = 0.0_dp

V_18_new = 0.0_dp 

V_28_new = 0.0_dp 

V_38_new = 0.0_dp 

!------------------------------------------------------------------
!(a) Compute recursively the alternative-specific value functions
!------------------------------------------------------------------
V_0_new = 0.0_dp


V_18_new = 0.0_dp/constv &
    + (1.0_dp - xi_21)*delta*( alpha_1*phi_grid(vi1) + beta_1*(1.0_dp - phi_grid(vi1)) )/constv &
	  *tau*( DMAX1(V_18_old(min_phi(vi1,1,2)),V_28_old(min_phi(vi1,1,2)),V_38_old(min_phi(vi1,1,2))) &
	  ! + 0.577215665_dp &
	  + DLOG(  DEXP(  ( (1.0_dp - lambda_11)*argu3(1, phi_grid(min_phi(vi1,1,2)), V_0_old(min_phi(vi1,1,2)) ) + lambda_11*value8_old(min_omega1low(1))     )/tau  &
	       - DMAX1(V_18_old(min_phi(vi1,1,2)),V_28_old(min_phi(vi1,1,2)),V_38_old(min_phi(vi1,1,2))) ) &
	  + DEXP( ( (1.0_dp - lambda_11)*argu3(1, phi_grid(min_phi(vi1,1,2)), V_18_old(min_phi(vi1,1,2)) ) + lambda_11*value8_old(min_omega1low(1))   )/tau  &
	       - DMAX1(V_18_old(min_phi(vi1,1,2)),V_28_old(min_phi(vi1,1,2)),V_38_old(min_phi(vi1,1,2)))  ) &
	  +DEXP(  (   (1.0_dp - lambda_11)*argu3(1, phi_grid(min_phi(vi1,1,2)), V_28_old(min_phi(vi1,1,2)) )  + lambda_11*value8_old(min_omega1low(1))       )/tau  &
	       - DMAX1(V_18_old(min_phi(vi1,1,2)),V_28_old(min_phi(vi1,1,2)),V_38_old(min_phi(vi1,1,2)))) &
	  +DEXP( ( (1.0_dp - lambda_11)*argu3(1, phi_grid(min_phi(vi1,1,2)), V_38_old(min_phi(vi1,1,2)) )  + lambda_11*value8_old(min_omega1low(1)) )/tau  &
	       - DMAX1(V_18_old(min_phi(vi1,1,2)),V_28_old(min_phi(vi1,1,2)),V_38_old(min_phi(vi1,1,2)))  ) )) &
	  + (1.0_dp - xi_21)*delta*(  (1.0_dp-alpha_1)*phi_grid(vi1) + (1.0_dp-beta_1)*(1.0_dp - phi_grid(vi1))  )/constv&
	  *tau*(DMAX1( V_18_old(min_phi(vi1,1,1)), V_28_old(min_phi(vi1,1,1)), V_38_old(min_phi(vi1,1,1))) &
		!+ 0.577215885_dp &
	   + DLOG( DEXP(   ( (1.0_dp - lambda_11)*argu3(1, phi_grid(min_phi(vi1,1,1)), V_0_old(min_phi(vi1,1,1)) ) + lambda_11*value8_old(min_omega1low(1))  ) /tau &
	       - DMAX1( V_18_old(min_phi(vi1,1,1)), V_28_old(min_phi(vi1,1,1)), V_38_old(min_phi(vi1,1,1)))  ) &
	  +DEXP( ( (1.0_dp - lambda_11)*argu3(1, phi_grid(min_phi(vi1,1,1)), V_18_old(min_phi(vi1,1,1)) ) + lambda_11*value8_old(min_omega1low(1))   )/tau  &
	      -  DMAX1( V_18_old(min_phi(vi1,1,1)), V_28_old(min_phi(vi1,1,1)), V_38_old(min_phi(vi1,1,1))) ) &
	  +DEXP( ( (1.0_dp - lambda_11)*argu3(1, phi_grid(min_phi(vi1,1,1)), V_28_old(min_phi(vi1,1,1)) )  + lambda_11*value8_old(min_omega1low(1))  )/tau &
	      -  DMAX1( V_18_old(min_phi(vi1,1,1)), V_28_old(min_phi(vi1,1,1)), V_38_old(min_phi(vi1,1,1))) ) &
	  +DEXP( (  (1.0_dp - lambda_11)*argu3(1, phi_grid(min_phi(vi1,1,1)), V_38_old(min_phi(vi1,1,1)) ) + lambda_11*value8_old(min_omega1low(1))  )/tau &
	      -  DMAX1( V_18_old(min_phi(vi1,1,1)), V_28_old(min_phi(vi1,1,1)), V_38_old(min_phi(vi1,1,1))) ))) 
	 




V_28_new =  (   (1.0_dp - delta)*(  d32 *(alpha_2*phi_grid(vi1)+beta_2*(1.0_dp-phi_grid(vi1)))   )  )/constv &
      + (1.0_dp - xi_25)*delta*( alpha_2*phi_grid(vi1) + beta_2*(1.0_dp - phi_grid(vi1)) )/constv &
	  *tau2*( DMAX1( V_18_old(min_phi(vi1,2,2)),V_28_old(min_phi(vi1,2,2)),V_38_old(min_phi(vi1,2,2))) &
		!+ 0.577215665_dp &
		+ DLOG(  DEXP( ((1.0_dp - lambda_12)*argu3(2, phi_grid(min_phi(vi1,2,2)), V_0_old(min_phi(vi1,2,2)) )  + lambda_12*value8_old(min_omega2low(1)) )/tau2 &
		   - DMAX1( V_18_old(min_phi(vi1,2,2)),V_28_old(min_phi(vi1,2,2)),V_38_old(min_phi(vi1,2,2))) ) &
		+ DEXP( ((1.0_dp - lambda_12)*argu3(2, phi_grid(min_phi(vi1,2,2)), V_18_old(min_phi(vi1,2,2)))    + lambda_12*value8_old(min_omega2low(1))    )/tau2 &
		   - DMAX1( V_18_old(min_phi(vi1,2,2)),V_28_old(min_phi(vi1,2,2)),V_38_old(min_phi(vi1,2,2)))) &
		+DEXP( ((1.0_dp - lambda_12)*argu3(2, phi_grid(min_phi(vi1,2,2)), V_28_old(min_phi(vi1,2,2))  )   + lambda_12*value8_old(min_omega2low(1))    )/tau2 &
		   - DMAX1( V_18_old(min_phi(vi1,2,2)),V_28_old(min_phi(vi1,2,2)),V_38_old(min_phi(vi1,2,2))) ) &
		+ DEXP( ((1.0_dp - lambda_12)*argu3(2, phi_grid(min_phi(vi1,2,2)), V_38_old(min_phi(vi1,2,2))  ) + lambda_12*value8_old(min_omega2low(1))    )/tau2 &
		   - DMAX1( V_18_old(min_phi(vi1,2,2)),V_28_old(min_phi(vi1,2,2)),V_38_old(min_phi(vi1,2,2))) ) ))  &
	  + (1.0_dp - xi_25)*delta*(  (1.0_dp-alpha_2)*phi_grid(vi1) + (1.0_dp-beta_2)*(1.0_dp - phi_grid(vi1))  )/constv &
	  *tau2*( DMAX1( V_18_old(min_phi(vi1,2,1)),V_28_old(min_phi(vi1,2,1)),V_38_old(min_phi(vi1,2,1))) &
		!+ 0.577215665_dp &
		+ DLOG(  DEXP(  ((1.0_dp - lambda_12)*argu3(2, phi_grid(min_phi(vi1,2,1)), V_0_old(min_phi(vi1,2,1)) ) + lambda_12*value8_old(min_omega2low(1))    )/tau2  &
		  -  DMAX1( V_18_old(min_phi(vi1,2,1)),V_28_old(min_phi(vi1,2,1)),V_38_old(min_phi(vi1,2,1))) ) &
		+ DEXP( ((1.0_dp - lambda_12)*argu3(2, phi_grid(min_phi(vi1,2,1)), V_18_old(min_phi(vi1,2,1)))  +  lambda_12*value8_old(min_omega2low(1))        )/tau2  &
		   - DMAX1( V_18_old(min_phi(vi1,2,1)),V_28_old(min_phi(vi1,2,1)),V_38_old(min_phi(vi1,2,1))) ) &
		+DEXP( ((1.0_dp - lambda_12)*argu3(2, phi_grid(min_phi(vi1,2,1)), V_28_old(min_phi(vi1,2,1))) + lambda_12*value8_old(min_omega2low(1))           )/tau2 &
		  -  DMAX1( V_18_old(min_phi(vi1,2,1)),V_28_old(min_phi(vi1,2,1)),V_38_old(min_phi(vi1,2,1)))  ) &
		+DEXP( ((1.0_dp - lambda_12)*argu3(2, phi_grid(min_phi(vi1,2,1)), V_38_old(min_phi(vi1,2,1))) +  lambda_12*value8_old(min_omega2low(1))         )/tau2 &
		  -  DMAX1( V_18_old(min_phi(vi1,2,1)),V_28_old(min_phi(vi1,2,1)),V_38_old(min_phi(vi1,2,1)))  ))) 
		





V_38_new =   (   (1.0_dp - delta)*(  d12 *(alpha_3*phi_grid(vi1)+beta_3*(1.0_dp-phi_grid(vi1)))  ) )/constv &
      + delta*( alpha_3*phi_grid(vi1) + beta_3*(1.0_dp - phi_grid(vi1)) )/constv &
	  *tau3*( DMAX1(V_18_old(min_phi(vi1,3,2)),V_28_old(min_phi(vi1,3,2)),V_38_old(min_phi(vi1,3,2)))  &
		!+ 0.577215665_dp &
		+ DLOG( DEXP( ((1.0_dp - lambda_13)*argu3(3, phi_grid(min_phi(vi1,3,2)), V_0_old(min_phi(vi1,3,2)) )  +  lambda_13*value8_old(min_omega3low(1))        )/tau3 &
		   - DMAX1(V_18_old(min_phi(vi1,3,2)),V_28_old(min_phi(vi1,3,2)),V_38_old(min_phi(vi1,3,2)))   ) &
		+ DEXP( ((1.0_dp - lambda_13)*argu3(3, phi_grid(min_phi(vi1,3,2)), V_18_old(min_phi(vi1,3,2)))   +  lambda_13*value8_old(min_omega3low(1))         )/tau3 &
		  -  DMAX1(V_18_old(min_phi(vi1,3,2)),V_28_old(min_phi(vi1,3,2)),V_38_old(min_phi(vi1,3,2)))  ) &
		+DEXP( ((1.0_dp - lambda_13)*argu3(3, phi_grid(min_phi(vi1,3,2)), V_28_old(min_phi(vi1,3,2)))  +  lambda_13*value8_old(min_omega3low(1))           )/tau3 &
		  -  DMAX1(V_18_old(min_phi(vi1,3,2)),V_28_old(min_phi(vi1,3,2)),V_38_old(min_phi(vi1,3,2)))   ) &
		+DEXP( ((1.0_dp - lambda_13)*argu3(3, phi_grid(min_phi(vi1,3,2)), V_38_old(min_phi(vi1,3,2)))  +  lambda_13*value8_old(min_omega3low(1))           )/tau3 &
		  -  DMAX1(V_18_old(min_phi(vi1,3,2)),V_28_old(min_phi(vi1,3,2)),V_38_old(min_phi(vi1,3,2)))   )  ))  &
	  + delta*(  (1.0_dp-alpha_3)*phi_grid(vi1) + (1.0_dp-beta_3)*(1.0_dp - phi_grid(vi1))  )/constv&
	  *tau3*( DMAX1(V_18_old(min_phi(vi1,3,1)),V_28_old(min_phi(vi1,3,1)),V_38_old(min_phi(vi1,3,1)))  &
		!+ 0.577215665_dp &
		+ DLOG(  DEXP( ((1.0_dp - lambda_13)*argu3(3, phi_grid(min_phi(vi1,3,1)), V_0_old(min_phi(vi1,3,1)) )    +  lambda_13*value8_old(min_omega3low(1))    )/tau3 &
		  -  DMAX1(V_18_old(min_phi(vi1,3,1)),V_28_old(min_phi(vi1,3,1)),V_38_old(min_phi(vi1,3,1)))  ) &
		+ DEXP( ((1.0_dp - lambda_13)*argu3(3, phi_grid(min_phi(vi1,3,1)), V_18_old(min_phi(vi1,3,1)))  +  lambda_13*value8_old(min_omega3low(1))         )/tau3 &
		 -  DMAX1(V_18_old(min_phi(vi1,3,1)),V_28_old(min_phi(vi1,3,1)),V_38_old(min_phi(vi1,3,1)))  ) &
         +DEXP( ((1.0_dp - lambda_13)*argu3(3, phi_grid(min_phi(vi1,3,1)), V_28_old(min_phi(vi1,3,1))) +  lambda_13*value8_old(min_omega3low(1))           )/tau3 &
		   - DMAX1(V_18_old(min_phi(vi1,3,1)),V_28_old(min_phi(vi1,3,1)),V_38_old(min_phi(vi1,3,1)))  ) &
	     +DEXP( ((1.0_dp - lambda_13)*argu3(3, phi_grid(min_phi(vi1,3,1)), V_38_old(min_phi(vi1,3,1))) + lambda_13*value8_old(min_omega3low(1))           )/tau3 &
		   - DMAX1(V_18_old(min_phi(vi1,3,1)),V_28_old(min_phi(vi1,3,1)),V_38_old(min_phi(vi1,3,1))) ))) 
	   


!-----------------------------------------------
!(b) Construct the Bellman operator as an array 
!-----------------------------------------------
B_operator(1) =  V_0_new
B_operator(2) =  V_18_new 
B_operator(3) =  V_28_new 
B_operator(4) =  V_38_new 


!------------------------------
!(c) Update the value function
!------------------------------
value8_new = MAXVAL(B_operator) 


!*********************************************************************
!(4) Update values, policy and alternative-specific value functions
!*********************************************************************
!IF (value_new >= PBAR) THEN
   value8(vi1) = value8_new   

    V_0(vi1)   = V_0_new    
    V_18(vi1)  = V_18_new  
    V_28(vi1)  = V_28_new  
    V_38(vi1)  = V_38_new 
	

END DO !end of the loop on the grid size



!***********************************************************************
!(5) Update value differences and iteration index (iterate on value old 
!    until value_old and value are at most equal to tol_level)
!***********************************************************************
diff = MAXVAL(DMAX1(V_18-V_18_old,V_28-V_28_old,V_38-V_38_old))

it = it+1

IF (diff<=tol_level .or. it>max_iter) THEN

PRINT*, 'exit: iteration', it
PRINT*, 'exit: convergence diff V cost', diff
PRINT*, ''
EXIT

END IF
  
END DO   !end of the loop on values 


!***************************************************************
!(B) Compute the value functions in the remaining tenure years
!***************************************************************
!---------------------------
!a. Undo the normalization
!---------------------------
V_0   = constv*V_0 + max_exp 
V_18   = constv*V_18 + max_exp 
V_28   = constv*V_28 + max_exp 
V_38   = constv*V_38 + max_exp 
value8 = constv*value8 + max_exp 

V_17 = 0.0_dp
V_27 = 0.0_dp
V_37 = 0.0_dp

V_16 = 0.0_dp
V_26 = 0.0_dp
V_36 = 0.0_dp


V_15 = 0.0_dp
V_25 = 0.0_dp
V_35 = 0.0_dp


V_14 = 0.0_dp
V_24 = 0.0_dp
V_34 = 0.0_dp


V_13 = 0.0_dp
V_23 = 0.0_dp
V_33 = 0.0_dp


V_12 = 0.0_dp
V_22 = 0.0_dp
V_32 = 0.0_dp


V_11 = 0.0_dp
V_21 = 0.0_dp
V_31 = 0.0_dp

value1 = 0.0_dp
value2 = 0.0_dp
value3 = 0.0_dp
value4 = 0.0_dp
value5 = 0.0_dp
value6 = 0.0_dp
value7 = 0.0_dp



!-------------
!b. Tenure 7
!-------------
PRINT*, 'Value function at tenure 7'



DO vi1 = 1, phi_grid_size

DO vl3 = 1, 3

DO vl4 = 1, 3

DO vl5 = 1, 3

DO vl6 = 1, 3



tenure3_L1=-1*(vl3==1)
tenure3_L2=-1*(vl3==2)
tenure3_L3=-1*(vl3==3)

tenure4_L1=-1*(vl4==1)
tenure4_L2=-1*(vl4==2)
tenure4_L3=-1*(vl4==3)

tenure5_L1=-1*(vl5==1)
tenure5_L2=-1*(vl5==2)
tenure5_L3=-1*(vl5==3)

tenure6_L1=-1*(vl6==1)
tenure6_L2=-1*(vl6==2)
tenure6_L3=-1*(vl6==3)



V_17(vi1,vl3,vl4,vl5,vl6) = (1.0_dp - delta)*(drift(1,7,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
    + (1.0_dp - xi_21)*delta*(alpha_1*phi_grid(vi1) + beta_1*(1.0_dp - phi_grid(vi1))) &
    * (  largeno    + DLOG( DEXP(V_0(min_phi(vi1,1,2))-largeno) + DEXP( V_18(min_phi(vi1,1,2))-largeno) + DEXP(V_28(min_phi(vi1,1,2))-largeno) + DEXP(V_38(min_phi(vi1,1,2))-largeno) )) &
    
    + (1.0_dp - xi_21)*delta*((1.0_dp-alpha_1)*phi_grid(vi1) + (1.0_dp-beta_1)*(1.0_dp - phi_grid(vi1)))&
	* (largeno    + DLOG( DEXP(V_0(min_phi(vi1,1,1))-largeno) + DEXP( V_18(min_phi(vi1,1,1))-largeno) + DEXP(V_28(min_phi(vi1,1,1))-largeno) + DEXP(V_38(min_phi(vi1,1,1))-largeno))) 
	  



V_27(vi1,vl3,vl4,vl5,vl6) = (1.0_dp - delta)*((  d23 +d32 )*(alpha_2*phi_grid(vi1)+beta_2*(1.0_dp-phi_grid(vi1))) &

    + drift(2,7,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_25)*delta*(alpha_2*phi_grid(vi1) + beta_2*(1.0_dp - phi_grid(vi1))) &
    * ( largeno    + DLOG( DEXP(V_0(min_phi(vi1,2,2))-largeno) + DEXP( V_18(min_phi(vi1,2,2))-largeno) + DEXP(V_28(min_phi(vi1,2,2))-largeno) + DEXP(V_38(min_phi(vi1,2,2))-largeno))) &
	
	+ (1.0_dp - xi_25)*delta*(  (1.0_dp-alpha_2)*phi_grid(vi1) + (1.0_dp-beta_2)*(1.0_dp - phi_grid(vi1))) &
	* ( largeno    + DLOG( DEXP(V_0(min_phi(vi1,2,1))-largeno) + DEXP( V_18(min_phi(vi1,2,1))-largeno) + DEXP(V_28(min_phi(vi1,2,1))-largeno) + DEXP(V_38(min_phi(vi1,2,1))-largeno))) 
	



V_37(vi1,vl3,vl4,vl5,vl6) = (1.0_dp - delta)*(d21 *(alpha_3*phi_grid(vi1)+beta_3*(1.0_dp-phi_grid(vi1))) &

    + drift(3,7,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_25 )*delta*(alpha_3*phi_grid(vi1) + beta_3*(1.0_dp - phi_grid(vi1))) &
    * ( largeno    + DLOG( DEXP(V_0(min_phi(vi1,3,2))-largeno) + DEXP( V_18(min_phi(vi1,3,2))-largeno) + DEXP(V_28(min_phi(vi1,3,2))-largeno) + DEXP(V_38(min_phi(vi1,3,2))-largeno))) &
	
	+ (1.0_dp - xi_25 )*delta*(  (1.0_dp-alpha_3)*phi_grid(vi1) + (1.0_dp-beta_3)*(1.0_dp - phi_grid(vi1))) &
	* ( largeno    + DLOG( DEXP(V_0(min_phi(vi1,3,1))-largeno) + DEXP( V_18(min_phi(vi1,3,1))-largeno) + DEXP(V_28(min_phi(vi1,3,1))-largeno) + DEXP(V_38(min_phi(vi1,3,1))-largeno))) 
	


END DO

END DO

END DO

END DO

END DO



!-------------
!c. Tenure 6
!-------------
PRINT*, 'Value function at tenure 6'



DO vi1 = 1, phi_grid_size

DO vl3 = 1, 3

DO vl4 = 1, 3

DO vl5 = 1, 3


tenure3_L1=-1*(vl3==1)
tenure3_L2=-1*(vl3==2)
tenure3_L3=-1*(vl3==3)

tenure4_L1=-1*(vl4==1)
tenure4_L2=-1*(vl4==2)
tenure4_L3=-1*(vl4==3)

tenure5_L1=-1*(vl5==1)
tenure5_L2=-1*(vl5==2)
tenure5_L3=-1*(vl5==3)



V_16(vi1,vl3,vl4,vl5) = (1.0_dp - delta)*( drift(1,6,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
    + (1.0_dp - xi_21)*delta*(alpha_1*phi_grid(vi1) + beta_1*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,1,2))-largeno) + DEXP( V_17(min_phi(vi1,1,2),vl3,vl4,vl5,1)-largeno) + DEXP(V_27(min_phi(vi1,1,2),vl3,vl4,vl5,1)-largeno) + DEXP(V_37(min_phi(vi1,1,2),vl3,vl4,vl5,1) -largeno))) &
    
    + (1.0_dp - xi_21)*delta*((1.0_dp-alpha_1)*phi_grid(vi1) + (1.0_dp-beta_1)*(1.0_dp - phi_grid(vi1)))&
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,1,1))-largeno) + DEXP( V_17(min_phi(vi1,1,1),vl3,vl4,vl5,1)-largeno) + DEXP(V_27(min_phi(vi1,1,1),vl3,vl4,vl5,1)-largeno) + DEXP(V_37(min_phi(vi1,1,1),vl3,vl4,vl5,1)-largeno))) 
	  



V_26(vi1,vl3,vl4,vl5) = (1.0_dp - delta)*(( d23 +d32 )*(alpha_2*phi_grid(vi1)+beta_2*(1.0_dp-phi_grid(vi1))) &

    + drift(2,6,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - cpen)*delta*(alpha_2*phi_grid(vi1) + beta_2*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,2,2))-largeno) + DEXP( V_17(min_phi(vi1,2,2),vl3,vl4,vl5,2)-largeno) + DEXP(V_27(min_phi(vi1,2,2),vl3,vl4,vl5,2)-largeno) + DEXP(V_37(min_phi(vi1,2,2),vl3,vl4,vl5,2)-largeno))) &
	
	+ (1.0_dp - cpen)*delta*(  (1.0_dp-alpha_2)*phi_grid(vi1) + (1.0_dp-beta_2)*(1.0_dp - phi_grid(vi1))) &
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,2,1))-largeno) + DEXP( V_17(min_phi(vi1,2,1),vl3,vl4,vl5,2)-largeno) + DEXP(V_27(min_phi(vi1,2,1),vl3,vl4,vl5,2)-largeno) + DEXP(V_37(min_phi(vi1,2,1),vl3,vl4,vl5,2)-largeno))) 
	
	
        
V_36(vi1,vl3,vl4,vl5) =  (1.0_dp - delta)*( DMAX1(d9,d28,d20)*(alpha_3*phi_grid(vi1)+beta_3*(1.0_dp-phi_grid(vi1))) &

    + drift(3,6,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - cpen )*delta*(alpha_3*phi_grid(vi1) + beta_3*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,3,2))-largeno) + DEXP( V_17(min_phi(vi1,3,2),vl3,vl4,vl5,3)-largeno) + DEXP(V_27(min_phi(vi1,3,2),vl3,vl4,vl5,3)-largeno) + DEXP(V_37(min_phi(vi1,3,2),vl3,vl4,vl5,3)-largeno))) &
	
	+ (1.0_dp - cpen )*delta*(  (1.0_dp-alpha_3)*phi_grid(vi1) + (1.0_dp-beta_3)*(1.0_dp - phi_grid(vi1))) &
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,3,1))-largeno) + DEXP( V_17(min_phi(vi1,3,1),vl3,vl4,vl5,3)-largeno) + DEXP(V_27(min_phi(vi1,3,1),vl3,vl4,vl5,3)-largeno) + DEXP(V_37(min_phi(vi1,3,1),vl3,vl4,vl5,3)-largeno))) 
	


END DO

END DO

END DO

END DO


!-------------
!d. Tenure 5
!-------------
PRINT*, 'Value function at tenure 5'



DO vi1 = 1, phi_grid_size

DO vl3 = 1, 3

DO vl4 = 1, 3



tenure3_L1=-1*(vl3==1)
tenure3_L2=-1*(vl3==2)
tenure3_L3=-1*(vl3==3)

tenure4_L1=-1*(vl4==1)
tenure4_L2=-1*(vl4==2)
tenure4_L3=-1*(vl4==3)



V_15(vi1,vl3,vl4) = (1.0_dp - delta)*(  drift(1,5,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
    + (1.0_dp - xi_21)*delta*(alpha_1*phi_grid(vi1) + beta_1*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,1,2))-largeno) + DEXP( V_16(min_phi(vi1,1,2),vl3,vl4,1)-largeno) + DEXP(V_26(min_phi(vi1,1,2),vl3,vl4,1)-largeno) + DEXP(V_36(min_phi(vi1,1,2),vl3,vl4,1)-largeno))) &
    
    + (1.0_dp - xi_21)*delta*((1.0_dp-alpha_1)*phi_grid(vi1) + (1.0_dp-beta_1)*(1.0_dp - phi_grid(vi1)))&
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,1,1))-largeno) + DEXP( V_16(min_phi(vi1,1,1),vl3,vl4,1)-largeno) + DEXP(V_26(min_phi(vi1,1,1),vl3,vl4,1)-largeno) + DEXP(V_36(min_phi(vi1,1,1),vl3,vl4,1)-largeno))) 
	  
	  


V_25(vi1,vl3,vl4) = (1.0_dp - delta)*((  d16 +d32 )*(alpha_2*phi_grid(vi1)+beta_2*(1.0_dp-phi_grid(vi1))) &

    + drift(2,5,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_34)*delta*(alpha_2*phi_grid(vi1) + beta_2*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,2,2))-largeno) + DEXP( V_16(min_phi(vi1,2,2),vl3,vl4,2)-largeno) + DEXP(V_26(min_phi(vi1,2,2),vl3,vl4,2)-largeno) + DEXP(V_36(min_phi(vi1,2,2),vl3,vl4,2)-largeno))) &
	
	+ (1.0_dp - xi_34)*delta*(  (1.0_dp-alpha_2)*phi_grid(vi1) + (1.0_dp-beta_2)*(1.0_dp - phi_grid(vi1))) &
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,2,1))-largeno) + DEXP( V_16(min_phi(vi1,2,1),vl3,vl4,2)-largeno) + DEXP(V_26(min_phi(vi1,2,1),vl3,vl4,2)-largeno) + DEXP(V_36(min_phi(vi1,2,1),vl3,vl4,2)-largeno))) 
	
	

V_35(vi1,vl3,vl4) =  (1.0_dp - delta)*(( DMAX1(d9,d28,d20) )*(alpha_3*phi_grid(vi1)+beta_3*(1.0_dp-phi_grid(vi1))) &

    + drift(3,5,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_34 )*delta*(alpha_3*phi_grid(vi1) + beta_3*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,3,2))-largeno) + DEXP( V_16(min_phi(vi1,3,2),vl3,vl4,3)-largeno) + DEXP(V_26(min_phi(vi1,3,2),vl3,vl4,3)-largeno) + DEXP(V_36(min_phi(vi1,3,2),vl3,vl4,3)-largeno))) &
	
	+ (1.0_dp - xi_34 )*delta*(  (1.0_dp-alpha_3)*phi_grid(vi1) + (1.0_dp-beta_3)*(1.0_dp - phi_grid(vi1))) &
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,3,1))-largeno) + DEXP( V_16(min_phi(vi1,3,1),vl3,vl4,3)-largeno) + DEXP(V_26(min_phi(vi1,3,1),vl3,vl4,3)-largeno) + DEXP(V_36(min_phi(vi1,3,1),vl3,vl4,3)-largeno))) 
	


END DO

END DO

END DO


!-------------
!e. Tenure 4
!-------------
PRINT*, 'Value function at tenure 4'


DO vi1 = 1, phi_grid_size

DO vl3 = 1, 3


tenure3_L1=-1*(vl3==1)
tenure3_L2=-1*(vl3==2)
tenure3_L3=-1*(vl3==3)


V_14(vi1,vl3) = (1.0_dp - delta)*( drift(1,4,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
    + (1.0_dp - xi_11)*delta*(alpha_1*phi_grid(vi1) + beta_1*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,1,2))-largeno) + DEXP( V_15(min_phi(vi1,1,2),vl3,1)-largeno) + DEXP(V_25(min_phi(vi1,1,2),vl3,1)-largeno) + DEXP(V_35(min_phi(vi1,1,2),vl3,1)-largeno))) &
    
    + (1.0_dp - xi_11)*delta*((1.0_dp-alpha_1)*phi_grid(vi1) + (1.0_dp-beta_1)*(1.0_dp - phi_grid(vi1)))&
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,1,1))-largeno) + DEXP( V_15(min_phi(vi1,1,1),vl3,1)-largeno) + DEXP(V_25(min_phi(vi1,1,1),vl3,1)-largeno) + DEXP(V_35(min_phi(vi1,1,1),vl3,1)-largeno))) 
	  
    


V_24(vi1,vl3) = (1.0_dp - delta)*(( d20 +d )*(alpha_2*phi_grid(vi1)+beta_2*(1.0_dp-phi_grid(vi1))) &

    + drift(2,4,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_14)*delta*(alpha_2*phi_grid(vi1) + beta_2*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,2,2))-largeno) + DEXP( V_15(min_phi(vi1,2,2),vl3,2)-largeno) + DEXP(V_25(min_phi(vi1,2,2),vl3,2)-largeno) + DEXP(V_35(min_phi(vi1,2,2),vl3,2)-largeno))) &
	
	+ (1.0_dp - xi_14)*delta*(  (1.0_dp-alpha_2)*phi_grid(vi1) + (1.0_dp-beta_2)*(1.0_dp - phi_grid(vi1))) &
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,2,1))-largeno) + DEXP( V_15(min_phi(vi1,2,1),vl3,2)-largeno) + DEXP(V_25(min_phi(vi1,2,1),vl3,2)-largeno) + DEXP(V_35(min_phi(vi1,2,1),vl3,2)-largeno))) 
	

   

V_34(vi1,vl3) =  (1.0_dp - delta)*( ( d27 + DMAX1(d9,d28, d20) )*(alpha_3*phi_grid(vi1)+beta_3*(1.0_dp-phi_grid(vi1))) &

    + drift(3,4,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_14 )*delta*(alpha_3*phi_grid(vi1) + beta_3*(1.0_dp - phi_grid(vi1))) &
    * (largeno    +DLOG( DEXP(V_0(min_phi(vi1,3,2))-largeno) + DEXP( V_15(min_phi(vi1,3,2),vl3,3)-largeno) + DEXP(V_25(min_phi(vi1,3,2),vl3,3)-largeno) + DEXP(V_35(min_phi(vi1,3,2),vl3,3)-largeno))) &
	
	+ (1.0_dp - xi_14 )*delta*(  (1.0_dp-alpha_3)*phi_grid(vi1) + (1.0_dp-beta_3)*(1.0_dp - phi_grid(vi1))) &
	* (largeno    +DLOG( DEXP(V_0(min_phi(vi1,3,1))-largeno) + DEXP( V_15(min_phi(vi1,3,1),vl3,3)-largeno) + DEXP(V_25(min_phi(vi1,3,1),vl3,3)-largeno) + DEXP(V_35(min_phi(vi1,3,1),vl3,3)-largeno))) 
	


END DO

END DO


!-------------
!f. Tenure 3
!-------------
PRINT*, 'Value function at tenure 3'


DO vi1 = 1, phi_grid_size



V_13(vi1) = (1.0_dp - delta)*( drift(1,3,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
    + (1.0_dp - xi_11 - plus)*delta*(alpha_1*phi_grid(vi1) + beta_1*(1.0_dp - phi_grid(vi1))) &
    * (largeno +   DLOG( DEXP(V_0(min_phi(vi1,1,2))-largeno) + DEXP( V_14(min_phi(vi1,1,2),1)-largeno) + DEXP(V_24(min_phi(vi1,1,2),1)-largeno) + DEXP(V_34(min_phi(vi1,1,2),1)-largeno))) &
    
    + (1.0_dp - xi_11 - plus)*delta*((1.0_dp-alpha_1)*phi_grid(vi1) + (1.0_dp-beta_1)*(1.0_dp - phi_grid(vi1)))&
	* (largeno +   DLOG( DEXP(V_0(min_phi(vi1,1,1))-largeno) + DEXP( V_14(min_phi(vi1,1,1),1)-largeno) + DEXP(V_24(min_phi(vi1,1,1),1)-largeno) + DEXP(V_34(min_phi(vi1,1,1),1)-largeno))) 
	  
  

V_23(vi1) = (1.0_dp - delta)*(( d28 +d )*(alpha_2*phi_grid(vi1)+beta_2*(1.0_dp-phi_grid(vi1))) &

    + drift(2,3,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_2 - plus)*delta*(alpha_2*phi_grid(vi1) + beta_2*(1.0_dp - phi_grid(vi1))) &
    * (largeno +DLOG( DEXP(V_0(min_phi(vi1,2,2))-largeno) + DEXP( V_14(min_phi(vi1,2,2),2)-largeno) + DEXP(V_24(min_phi(vi1,2,2),2)-largeno) + DEXP(V_34(min_phi(vi1,2,2),2)-largeno))) &
	
	+ (1.0_dp - xi_2 - plus)*delta*(  (1.0_dp-alpha_2)*phi_grid(vi1) + (1.0_dp-beta_2)*(1.0_dp - phi_grid(vi1))) &
	* (largeno +DLOG( DEXP(V_0(min_phi(vi1,2,1))-largeno) + DEXP( V_14(min_phi(vi1,2,1),2)-largeno) + DEXP(V_24(min_phi(vi1,2,1),2)-largeno) + DEXP(V_34(min_phi(vi1,2,1),2)-largeno))) 

     

V_33(vi1) =  (1.0_dp - delta)*((-d27)*(alpha_3*phi_grid(vi1)+beta_3*(1.0_dp-phi_grid(vi1))) &

    + drift(3,3,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_3)*delta*(alpha_3*phi_grid(vi1) + beta_3*(1.0_dp - phi_grid(vi1))) &
    * (largeno +DLOG( DEXP(V_0(min_phi(vi1,3,2))-largeno) + DEXP( V_14(min_phi(vi1,3,2),3)-largeno) + DEXP(V_24(min_phi(vi1,3,2),3)-largeno) + DEXP(V_34(min_phi(vi1,3,2),3)-largeno))) &
	
	+ (1.0_dp - xi_3)*delta*(  (1.0_dp-alpha_3)*phi_grid(vi1) + (1.0_dp-beta_3)*(1.0_dp - phi_grid(vi1))) &
	* (largeno +DLOG( DEXP(V_0(min_phi(vi1,3,1))-largeno) + DEXP( V_14(min_phi(vi1,3,1),3)-largeno) + DEXP(V_24(min_phi(vi1,3,1),3)-largeno) + DEXP(V_34(min_phi(vi1,3,1),3)-largeno))) 
	

END DO



!-------------
!g. Tenure 2
!-------------
PRINT*, 'Value function at tenure 2'


DO vi1 = 1, phi_grid_size


V_12(vi1) = (1.0_dp - delta)*(d10 *(alpha_1*phi_grid(vi1)+beta_1*(1.0_dp-phi_grid(vi1))) & 
       
    + drift(1,2,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
    + (1.0_dp - xi_1)*delta*(alpha_1*phi_grid(vi1) + beta_1*(1.0_dp - phi_grid(vi1))) &
    * (largeno +DLOG( DEXP(V_0(min_phi(vi1,1,2))-largeno) + DEXP( V_13(min_phi(vi1,1,2))-largeno) + DEXP(V_23(min_phi(vi1,1,2))-largeno) + DEXP(V_33(min_phi(vi1,1,2))-largeno))) &
    
    + (1.0_dp - xi_1)*delta*((1.0_dp-alpha_1)*phi_grid(vi1) + (1.0_dp-beta_1)*(1.0_dp - phi_grid(vi1)))&
	* (largeno +DLOG( DEXP(V_0(min_phi(vi1,1,1))-largeno) + DEXP( V_13(min_phi(vi1,1,1))-largeno) + DEXP(V_23(min_phi(vi1,1,1))-largeno) + DEXP(V_33(min_phi(vi1,1,1))-largeno))) 
	  



V_22(vi1) = (1.0_dp - delta)*((d9 +d )*(alpha_2*phi_grid(vi1)+beta_2*(1.0_dp-phi_grid(vi1))) &

    + drift(2,2,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_2)*delta*(alpha_2*phi_grid(vi1) + beta_2*(1.0_dp - phi_grid(vi1))) &
    * (largeno +DLOG( DEXP(V_0(min_phi(vi1,2,2))-largeno) + DEXP( V_13(min_phi(vi1,2,2))-largeno) + DEXP(V_23(min_phi(vi1,2,2))-largeno) + DEXP(V_33(min_phi(vi1,2,2))-largeno))) &
	
	+ (1.0_dp - xi_2)*delta*(  (1.0_dp-alpha_2)*phi_grid(vi1) + (1.0_dp-beta_2)*(1.0_dp - phi_grid(vi1))) &
	* (largeno +DLOG( DEXP(V_0(min_phi(vi1,2,1))-largeno) + DEXP( V_13(min_phi(vi1,2,1))-largeno) + DEXP(V_23(min_phi(vi1,2,1))-largeno) + DEXP(V_33(min_phi(vi1,2,1))-largeno))) 


     

V_32(vi1) =  (1.0_dp - delta)*((-d27) *(alpha_3*phi_grid(vi1)+beta_3*(1.0_dp-phi_grid(vi1))) &

    + drift(3,2,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_3)*delta*(alpha_3*phi_grid(vi1) + beta_3*(1.0_dp - phi_grid(vi1))) &
    * (largeno +DLOG( DEXP(V_0(min_phi(vi1,3,2))-largeno) + DEXP( V_13(min_phi(vi1,3,2))-largeno) + DEXP(V_23(min_phi(vi1,3,2))-largeno) + DEXP(V_33(min_phi(vi1,3,2))-largeno))) &
	
	+ (1.0_dp - xi_3)*delta*(  (1.0_dp-alpha_3)*phi_grid(vi1) + (1.0_dp-beta_3)*(1.0_dp - phi_grid(vi1))) &
	* (largeno +DLOG( DEXP(V_0(min_phi(vi1,3,1))-largeno) + DEXP( V_13(min_phi(vi1,3,1))-largeno) + DEXP(V_23(min_phi(vi1,3,1))-largeno) + DEXP(V_33(min_phi(vi1,3,1))-largeno))) 
	
	


END DO



!-------------
!h. Tenure 1
!-------------
PRINT*, 'Value function at tenure 1'


DO vi1 = 1, phi_grid_size


V_11(vi1) = (1.0_dp - delta)*( drift(1,1,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
    + (1.0_dp - xi_1)*delta*(alpha_1*phi_grid(vi1) + beta_1*(1.0_dp - phi_grid(vi1))) &
    * (largeno +DLOG( DEXP(V_0(min_phi(vi1,1,2))-largeno) + DEXP( V_12(min_phi(vi1,1,2))-largeno) + DEXP(V_22(min_phi(vi1,1,2))-largeno) + DEXP(V_32(min_phi(vi1,1,2))-largeno))) &
    
    + (1.0_dp - xi_1)*delta*((1.0_dp-alpha_1)*phi_grid(vi1) + (1.0_dp-beta_1)*(1.0_dp - phi_grid(vi1)))&
	* (largeno +DLOG( DEXP(V_0(min_phi(vi1,1,1))-largeno) + DEXP( V_12(min_phi(vi1,1,1))-largeno) + DEXP(V_22(min_phi(vi1,1,1))-largeno) + DEXP(V_32(min_phi(vi1,1,1))-largeno))) 
	  
	  

V_21(vi1) = (1.0_dp - delta)*( drift(2,1,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_2)*delta*(alpha_2*phi_grid(vi1) + beta_2*(1.0_dp - phi_grid(vi1))) &
    * (largeno +DLOG( DEXP(V_0(min_phi(vi1,2,2))-largeno) + DEXP( V_12(min_phi(vi1,2,2))-largeno) + DEXP(V_22(min_phi(vi1,2,2))-largeno) + DEXP(V_32(min_phi(vi1,2,2))-largeno))) &
	
	+ (1.0_dp - xi_2)*delta*(  (1.0_dp-alpha_2)*phi_grid(vi1) + (1.0_dp-beta_2)*(1.0_dp - phi_grid(vi1))) &
	* (largeno +DLOG( DEXP(V_0(min_phi(vi1,2,1))-largeno) + DEXP( V_12(min_phi(vi1,2,1))-largeno) + DEXP(V_22(min_phi(vi1,2,1))-largeno) + DEXP(V_32(min_phi(vi1,2,1))-largeno))) 



V_31(vi1) =  (1.0_dp - delta)*((-d27)*(alpha_3*phi_grid(vi1)+beta_3*(1.0_dp-phi_grid(vi1))) &

    + drift(3,1,1,1,tenure3_L1,tenure3_L2,tenure3_L3,tenure4_L1,tenure4_L2,tenure4_L3,tenure5_L1,tenure5_L2,tenure5_L3,tenure6_L1,tenure6_L2,tenure6_L3))  &
                       
    + (1.0_dp - xi_3)*delta*(alpha_3*phi_grid(vi1) + beta_3*(1.0_dp - phi_grid(vi1))) &
    * (largeno +DLOG( DEXP(V_0(min_phi(vi1,3,2))-largeno) + DEXP( V_12(min_phi(vi1,3,2))-largeno) + DEXP(V_22(min_phi(vi1,3,2))-largeno) + DEXP(V_32(min_phi(vi1,3,2))-largeno))) &
	
	+ (1.0_dp - xi_3)*delta*(  (1.0_dp-alpha_3)*phi_grid(vi1) + (1.0_dp-beta_3)*(1.0_dp - phi_grid(vi1))) &
	* (largeno +DLOG( DEXP(V_0(min_phi(vi1,3,1))-largeno) + DEXP( V_12(min_phi(vi1,3,1))-largeno) + DEXP(V_22(min_phi(vi1,3,1))-largeno) + DEXP(V_32(min_phi(vi1,3,1))-largeno))) 
	

END DO


phi_1_star = 0.5_dp 
phi_2_star = 0.5_dp 
phi_3_star = 0.5_dp

phi_1_grid_point = 10
phi_1_grid_point_max = 20


END SUBROUTINE value_infinite



!*************************************************************************
!(C) Function to make the value function continuous used in prob_wage
!*************************************************************************
FUNCTION value_cont(vfct,phi)

REAL(dp), DIMENSION(phi_grid_size,totsample), INTENT(IN):: vfct
REAL(dp), INTENT(IN):: phi
REAL(dp):: value_cont
INTEGER:: index


!-------------------------------------------------------------------
!(a) No need for approximation if the belief value is on the grid
!-------------------------------------------------------------------
DO index = 1, phi_grid_size 
 IF (phi==phi_grid(index)) THEN
  value_cont = vfct(index,1)
 END IF
END DO


!------------------------------------------------------------------------
!(b) Linear interpolation if the belief value falls btw two consecutive
!    values on the grid: y=y1+[(y2-y1)/(x2-x1)](x-x1) 
!------------------------------------------------------------------------
DO index = 2, phi_grid_size  
 IF (phi_grid(index-1) < phi .and. phi < phi_grid(index)) THEN

   value_cont = vfct(index-1,1)+((vfct(index,1)-vfct(index-1,1))/(phi_grid(index)-phi_grid(index-1)))&
              *(phi-phi_grid(index-1))

 END IF
END DO

               
END FUNCTION value_cont


END MODULE mod_values
